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25  Abstract 

26  The  non-hydrostatic  (NH)  compressible  Euler  equations  of  dry  atmosphere  are  solved  in 

27  a  simplified  two  dimensional  (2D)  slice  (X-Z)  framework  employing  a  spectral  element 

28  method  (SEM)  for  the  horizontal  discretization  and  a  finite  difference  method  (EDM)  for  the 

29  vertical  discretization.  The  SEM  uses  high-order  nodal  basis  functions  associated  with 

30  Eagrange  polynomials  based  on  Gauss-Eobatto-Legendre  (GEE)  quadrature  points.  The  EDM 

31  employs  a  third-order  upwind  biased  scheme  for  the  vertical  flux  terms  and  a  centered  finite 

32  difference  scheme  for  the  vertical  derivative  terms  and  quadrature.  The  Euler  equations  used 

33  here  are  in  a  flux  form  based  on  the  hydrostatic  pressure  vertical  coordinate,  which  are  the 

34  same  as  those  used  in  the  Weather  Research  and  Eorecasting  (WRE)  model,  but  a  hybrid 

35  sigma-pressure  vertical  coordinate  is  implemented  in  this  model.  We  verified  the  model  by 

36  conducting  widely  used  standard  benchmark  tests:  the  inertia-gravity  wave,  rising  thermal 

37  bubble,  density  current  wave,  and  linear  hydrostatic  mountain  wave.  The  numerical  results 

38  demonstrate  that  the  horizontally  spectral  element  vertically  finite  difference  model  is 

39  accurate  and  robust.  By  using  the  2D  slice  model,  we  effectively  show  that  the  combined 

40  spatial  discretization  method  of  the  spectral  element  and  finite  difference  method  in  the 

41  horizontal  and  vertical  directions,  respectively,  offers  a  viable  method  for  the  development  of 

42  a  NH  dynamical  core.  The  present  core  provides  a  practical  framework  for  further 

43  development  of  three-dimensional  (3D)  non-hydrostatic  compressible  atmospheric  models. 

44 
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1.  Introduction 


46  There  is  a  growing  interest  in  developing  highly  sealable  dynamieal  eores  using 

47  numerieal  algorithms  under  petaseale  eomputers  with  many  eores  (with  the  goal  of  exaseale 

48  eomputing  just  around  the  eorner).  The  speetral  element  method  (SEM)  has  been  known  as 

49  one  of  the  most  promising  methods  with  high  effieieney  and  aeouraey  (Taylor  et  al.  1997; 

50  Giraldo  2001;  Thomas  and  Loft  2002).  SEM  is  loeal  in  nature  beeause  of  having  a  large  on- 

51  proeessor  operation  eount  (Kelly  and  Graldo,  2012).  The  SEM  aehieves  this  high  level  of 

52  sealability  by  deeomposing  the  physieal  domain  into  smaller  pieees  with  a  small 

53  eommunieation  steneil.  Also  SEM  has  been  shown  to  be  very  attraetive  in  aehieving  high- 

54  order  aeeuraey  and  geometrieal  flexibility  on  the  sphere  (Taylor  et  al.  1997;  Giraldo  2001; 

55  Giraldo  et  al.  2004). 

56  To  date,  the  SEM  has  been  sueeessfully  implemented  in  atmospherie  modeling  sueh  as 

57  in  the  Community  Atmosphere  Model  -  speetral  element  dynamieal  eore  (CAM-SE) 

58  (Thomas  and  Loft  2005)  and  the  sealable  speetral  element  Eulerian  atmospherie  model  (SEE- 

59  AM)  (Giraldo  and  Rosmond,  2004).  These  models  eonsider  the  primitive  hydrostatie 

60  equations  on  global  grid  meshes  sueh  as  a  eubed-sphere  tiled  with  quadrilateral  elements 

61  using  SEM  in  the  horizontal  diseretization  and  the  finite  differenee  method  (EDM)  in  the 

62  vertieal.  The  robustness  of  the  SEM  has  been  illustrated  through  three-dimensional  dry 

63  dynamieal  test  eases  (Thomas  and  Loft  2005;  Giraldo  and  Rosmond  2004;  Giraldo  2005; 

64  Taylor  et  al.  2007;  Lauritzen  et  al.  2010). 

65  The  ultimate  objeetive  of  our  study  is  to  build  a  3D  non-hydrostatie  (NH)  model  based 

66  on  the  eompressible  Navier-Stokes  equations  using  the  eombined  horizontally  SEM  and 

67  vertieally  EDM.  Sinee  testing  a  3D  NH  model  requires  a  huge  amount  of  eomputing 

68  resourees,  studying  the  feasibility  of  our  approaeh  in  2D  is  an  attraetive  alternatively  to  the 
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69  development  of  a  fully  3D  model.  This  is  the  ease  because  a  2D  slice  model  effectively  can 

70  test  the  practical  issues  resulting  from  the  vertical  discretization  and  time  integration,  prior  to 

71  the  construction  of  a  full  3D  model.  Although  we  could  also  discretize  the  vertical  direction 

72  with  SEM  (as  is  proposed  in  Kelly  and  Giraldo  2012  and  Giraldo  et  al.  2013),  we  choose  to 

73  use  a  conservative  flux-form  finite-difference  method  for  discretization  in  the  vertical 

74  direction,  which  provides  an  easy  way  for  coupling  the  dynamics  and  existing  physics 

75  packages. 

76  We  have  developed  a  dry  2D  NH  compressible  Euler  model  based  on  SEM  along  the  x- 

77  direction  and  EDM  along  the  z-direction  for  this  purpose.  Hereafter,  this  is  simply  referred  to 

78  as  the  2DNH  model.  We  adopt  the  governing  equation  formulation  proposed  by  Skamarock 

79  and  Klemp  (2008)  (hereafter,  SK08)  which  is  used  in  the  Weather  Research  and  Eorecasting 

80  (WRE)  Model.  The  Euler  equations  are  in  flux  form  based  on  the  hydrostatic  pressure  vertical 

81  coordinate.  In  SK08,  the  terrain-following  sigma-pressure  coordinate  is  used,  but  here  we 

82  employ  a  hybrid  sigma-pressure  vertical  coordinate.  Park  et  al.  (2013)  (hereafter,  PK13) 

83  provides  a  clue  for  the  equation  set  in  the  hybrid  sigma-pressure  in  their  appendix,  in  which 

84  the  hybrid  sigma-pressure  coordinate  is  applied  to  the  hydrostatic  primitive  equations  and  can 

85  be  modified  exactly  to  the  sigma-pressure  coordinate  at  the  level  of  the  actual  coding 

86  implementation.  Also,  we  built  the  2DNH  model  using  a  time-split  third-order  Runge-Kutta 

87  (RK3)  for  the  time  discretization,  which  has  been  shown  to  work  effectively  in  the  WRE 

88  model.  We  keep  the  temporal  discretization  of  the  model  as  similar  as  possible  to  the  WRE 

89  model  in  order  to  more  directly  the  discern  the  differences  related  to  the  discrete  spatial 

90  operators  between  the  two  models.  This  provides  robust  tools  for  development  and 

91  verification  of  the  2DNH  model. 

92  In  this  paper,  we  show  the  feasibility  of  the  2DNH  model  by  conducting  conventional 

93  benchmark  test  cases  as  well  as  focusing  on  the  description  of  the  numerical  scheme  for  the 
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94  spatial  discretization.  We  verify  the  2DNH  by  analyzing  four  test  cases:  the  inertia-gravity 

95  wave,  rising  thermal  bubble,  density  eurrent  wave,  and  linear  hydrostatie  mountain  wave. 

96  The  organization  of  this  paper  is  as  follows.  The  next  seetion  deseribes  the  governing 

97  equations  with  a  definition  of  the  prognostie  and  diagnostic  variables  used  in  our  model,  in 

98  which  we  present  essential  ehanges  from  SK08.  Section  3  contains  the  description  of  the 

99  temporal  and  spatial  diseretization  ineluding  the  speetral  element  formulation.  In  Sec.  4,  we 

100  present  the  results  of  the  2DNH  model  using  all  four  test  cases.  Finally  in  See.  5  we 

101  summarize  the  paper  and  propose  future  direetions. 

102 

103  2.  Governing  equations 

104  We  adopt  the  governing  equation  formulation  of  SK08.  Here  we  implement  the  hybrid 

105  sigma-pressure  eoordinate  reported  by  PK13  whieh  eonsidered  only  the  hydrostatie  primitive 

106  equation.  The  hybrid  sigma  pressure  eoordinate  is  defined  with  rj  e  |^0, 1^  as 

107  p^=B{vi)[p^-p)  +  [ri-B{ri)][p^-p)  +  p^  (1) 

108  where  is  the  hydrostatie  pressure  of  dry  air,  Bij])  is  the  relative  weighting  of  the 

109  terrain-following  coordinate  versus  the  normalized  pressure  coordinate,  p^ ,  p^ ,  and  p^are 

110  the  hydrostatic  surface  pressure  of  dry  air,  the  top  level  pressure,  and  a  reference  sea  level 

1 1 1  pressure,  respeetively.  A  more  detailed  deseription  of  the  hybrid  sigma  pressure  eoordinate 

1 12  ean  be  found  in  the  Appendix  of  PK13.  The  definition  of  the  flux  variables  are 

113  {y^,\NM,Q)  =  (2) 

114  where  and  w  are  the  veloeities  in  the  horizontal  and  vertical  directions, 

dll 

115  respectively,  77  =  —  is  the  77 -eoordinate  (eontravariant)  vertical  velocity,  0  is  the 

dt 
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117 

118 

119 


120 


121 


122 


123 


potential  temperature,  and  is  the  mass  of  the  dry  air  in  the  layers  defined  as 


1- 


3e(?7) 

dr] 


(Po-Pt)' 


The  flux-form  Euler  equations  for  dry  atmosphere  are  expressed  as 


dV, 

dt 


=  -Pd 


—  -P' 

dv 


drj 


m 

dt 


=  g 


dl] 


Pa 


dt  dt 


dp 


\  dP  J 


dB{i])  dp' 


dQ 


- ^  =  -V  -1/  -  — . 

dr]  dt  ”  ^  dr] 


d<p' 

dt 


Pd 


30 

dt 


d(QO) 

dr] 


(3) 

(4) 

(5) 

(6) 

(V) 

(8) 


124  where  0  is  the  geopotential,  is  the  inverse  density  for  dry  air,  and  F  r  and 

125  represent  foreing  terms  of  the  Coriolis  and  curvature  which  we  ignore  for  simplicity.  In  Eqs. 

126  (4)-(8),  the  governing  equations  are  described  with  perturbation  variables  such  as 

127  p  =  p{z)  +  p',  <p  =  0{z)  +  <p',  =  ajz~)  +  a'^  ,  and  p^  =  p^{x,y)  +  p'^  where  the 

128  variables  denoted  by  overbars  are  reference  state  variables  that  satisfy  hydrostatic  balance. 

129  Eor  completeness,  the  diagnostic  relation  for  Q  is  given  by  integrating  Eq.  (6) 

130  vertically  from  the  surface  (r]  =  1)  to  the  material  surface  as 


131 


dB(r])  dp' 
dr]  dt 


+  V  -1/^77, 


(9) 


132 


where 


dt 


is  obtained  by  integrating  vertically  Eq.  (6)  for  the  surface  (r]  =  1)  to  the  top 
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133  (rj  =  0)  using  a  no-flux  boundary  condition  such  as  Q|  ^  ^  =  0  and  the  specification  of 

134  the  vertical  coordinate  sueh  as  Birj  =  1)  =  1  and  Birj  =  0)  =  0  as 


135 

136 

137 

138 

139 


(10) 


The  above  equation  allows  /O'  to  be  evolved  forward  in  time  where  we  then  compute 
directly  from  Eq.  (5).  The  diagnostic  relation  for  the  dry  inverse  density  is  given  as 


dri 


^dl^i 


d^d  ’ 


(11) 


and  the  full  pressure  for  dry  atmosphere  is 


140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 


P  =Po 


(12) 


This  concludes  the  deseription  of  the  governing  equations  used  in  our  model;  in  the  next 
section  we  deseribe  the  discretization  of  the  continuous  form  of  the  governing  equations  that 
are  used  in  our  model. 


3.  Discretization 


a.  Spatial  diseretization 
1)  Horizontal  direetion 

For  a  given  rj  level,  we  discretize  the  horizontal  operators  using  the  SEM.  Therefore  in 
2D  (X-Z)  sliee  framework  we  foeus  on  the  SEM  disere te  gradient  operator  for  ID  (x).  In 

SEM,  we  approximate  the  solution  in  non-overlapping  elements  as 

w+i 

=  (13) 

k=\ 

where  x^,  represents  A/ -1-1  grid  points  that  correspond  to  the  Gauss-Eobatto-Eegendre 
(GEE)  points  and  V^^(x)  are  the  A/ th-order  Eagrange  polynomials  based  on  the  GEE  points. 
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153 

154 

155 

156 

157 


158 


159 

160 


161 


162 

163 

164 

165 

166 

167 


168 

169 

170 

171 


It  is  noteworthy  that  the  i//"^  have  the  cardinal  property,  i.e.,  they  can  be  represented  as 
Kronecker  delta  functions  where  i//"^  are  zero  at  all  nodal  points  except  (but  are 
allowed  to  oscillate  between  nodal  points). 

The  GLL  points  in  a  reference  coordinate  system  ^g[-1,+1]  and  the  associated 
quadrature  weights  , 

COi^J  =  - - r 

/V(/V  +  1) 

are  introduced  for  the  Gaussian  quadrature; 

I  ?  dir  = 

/=0 

dx 

where  are  the  N  th-order  Legendre  polynomials,  J  =  —  is  the  transformation 

/v  ^  ff 


PA) 


(14) 


Jacobian,  and  Q*"  represents  the  non-overlapping  elements. 

We  now  introduce  the  polynomial  expansions  into  our  governing  equations  in  the  form 


of 


dt 


(16) 


multiply  by  the  basis  function  as  a  test  function,  and  integrate  to  yield  a  system  of  ordinary 
differential  equations  as  such 


A/+1 

n=] 


dQ, 


e 

nk 


dt 


(  A/+1 

T  p 


Vn=1 


k  d^.. 


(17) 


where  k  =  1 2,L  ,N  +1,  is  the  element  based  mass  matrix  given  as 


^"nk=\^  ¥n¥k  d^  =  CO^\j,\S^,, 


(18) 


and  the  right-hand  side  of  Eq.  (18)  is  evaluated  using  Gaussian  quadrature  of  Eq.  (16).  It  is 
noted  that  using  GEE  points  for  both  interpolation  and  integration  results  in  a  diagonal  mass 
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172 


matrix  M^,  ,  which  means  that  the  inversion  of  the  mass  matrix  is  trivial. 

The  horizontal  derivatives  included  in  the  right-hand  side  of  Eq.  (17)  are  evaluated  using 
the  analytie  derivatives  of  the  basis  funetions  as  follows 


173 

174 

175 


176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 


187 


188 

189 

190 


191 


ax  ~  d^dX~ 


/<=! 


ax 


(19) 


Note  that  the  non-differential  operations,  such  as  cross  products,  are  eomputed  direetly  at  grid 
points  since  we  use  nodal  basis  functions  associated  with  Lagrange  polynomials  based  on  the 
GEL  points.  In  order  to  satisfy  the  equations  globally,  we  use  the  direct  stiffness  summation 
(DSS)  operation.  Eor  a  more  detailed  deseription  of  the  SEM,  see  Giraldo  and  Rosmond 
(2004),  Giraldo  and  Restelli  (2008),  and  Kelly  and  Giraldo  (2012). 


2)  Vertieal  direetion 

Using  a  Lorenz  staggering,  the  variables  1/^ ,  0 ,  ju ,  a ,  p  are  at  layer  midpoints 
denoted  by  k  -'X2,K  ,K  where  K  is  the  total  number  of  layers,  and  the  variables  W , 

(j)  live  at  layer  interfaees  defined  by  k  +  ^,  k  —0,XK  ,K  ,  so  that  and 

^V2  “  ^Bottom  ~  ^  deseribes  the  grid  points  and  the  alloeation  of  the  variables.  Here,  we 

evaluate  the  vertieal  adveetion  terms  (  — - -  ,  — - -  ,  and  — - -  )  and  vertieal 

dri  dri  dri 


derivative  terms  ( ,  and  — ).  The  former  is  discretized  using  the  third-order  upwind 

dri  dri 

biased  diseretization  in  Hundsdorfer  et  al.  (1995)  whieh  is  given  as 


dri 


4_2  84-1  +  84^1  Lz,  ^  3ign(^)  ,  (20) 


1 2A77 


1 2Ari 


where  f  corresponds  to  the  flux  sueh  as  ,  and  A77  =  ri^  -  rii^_y„  is  the  thickness  of 
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192  the  layer.  The  latter  is  discretized  by  the  centered  finite  difference.  Likewise  the  vertical 

193  discretization  quadrature  rules  for  the  calculations  of  Eqs.  (9)  and  (10)  follow  the  finite 

194  difference  discretization  naturally. 

195 

196  b.  Temporal  discretization 

197  For  integrating  the  equations,  we  use  the  time-split  RK3  integration  technique  following 

198  the  strategy  of  SK08,  in  which  low-frequency  modes  due  to  advective  forcings  are  explicitly 

199  advanced  using  a  large  time  step  of  the  RK3  scheme,  but  high-frequency  modes  are 

200  integrated  over  smaller  time  steps  using  an  explicit  forward-backward  time  integration 

201  scheme  for  the  horizontally  propagating  acoustic/gravity  waves  and  a  fully  implicit  scheme 

202  for  vertically  propagating  acoustic  waves  and  buoyancy  oscillations  (Klemp  et  al.  2007)  . 

203  This  technique  has  been  shown  to  work  effectively  within  numerous  nonhydrostatic  models 

204  including  the  WRF  model  (Skamarock  et  al.  2008),  the  Model  for  Prediction  Across  Scales 

205  (MPAS)  (Skamarock  et  al.  2012),  and  the  Nonhydrostatic  Icosahedral  Atmospheric  Model 

206  (NICAM)  (Satoh  et  al.  2008). 

207  It  is  noted  that  in  the  procedure  of  the  time-split  RK3  integration,  the  difference  between 

208  the  approach  used  in  this  paper  and  SK08  comes  from  the  vertical  coordinate.  Since  we  use 

209  the  hybrid  sigma-pressure  coordinate,  the  equation  for  p'^  (Eq.  (6))  should  be  first  stepped 

210  forward  in  time  using  the  forward-backward  differencing  on  the  small  time  steps,  then 

211  can  be  computed  directly  from  the  specification  of  the  vertical  coordinate  in  Eq.  (9)  and 

212  is  obtained  from  the  vertical  integration. 

213 

214  4.  Test  cases 

215  We  validate  the  2DNH  model  on  four  test  cases  of  the  linear  hydrostatic  mountain  wave. 
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216  density  current,  inertia-gravity  wave,  and  rising  thermal  bubble  experiments.  All  cases  but  the 

217  mountain  wave  experiment  do  not  have  analytic  solutions.  Therefore,  for  the  mountain  wave 

218  experiment,  numerical  results  of  the  2DNH  model  are  compared  to  analytic  solutions  (Durran 

219  and  Klemp  1983),  and  for  the  other  experiments,  we  compare  our  results  to  the  results  of 

220  other  published  papers. 

221  It  should  be  mentioned  that  the  horizontal  SEM  formulation  is  able  to  utilize  arbitrary 

222  order  polynomials  per  element  to  represent  the  discrete  spatial  operators,  but  in  this  paper  all 

223  the  results  presented  use  either  5th  or  8th  order  polynomials.  The  averaged  horizontal  grid 

224  spacing  is  defined  as 


225 


Ax  = 


n=^ 


N 


(21) 


226  where  Ax^  is  the  internal  grid  spacing  within  the  element  which  is  regularly  spaced  in  the 

227  domain  and  N  is  the  number  of  the  interval  associated  with  irregularly  spaced  GEL 

228  quadrature  points  which  is  equivalent  to  the  order  of  the  basis  polynomials.  The  average 

229  vertical  grid  spacing  is  defined  in  the  same  way  of  Eq.  (24).  Below,  we  use  this  convention  to 

230  define  the  grid  resolution. 

231 

232  a.  Einear  hydrostatic  mountain  wave  test 

233  In  order  to  verify  the  2DNH’s  feasibility  to  treat  surface  elevations  associated  with  the 

234  vertical  terrain-following  coordinate,  we  simulate  the  linear  hydrostatic  mountain  wave  test 

235  introduced  by  Durran  and  Klemp  (1983)  (hereafter,  DK83)  in  which  the  analytic  steady-state 

236  solution  is  provided  by  using  a  single-peaked  mountain  with  uniform  zonal  wind.  To  compare 

237  our  results  with  the  analytic  and  numerical  solution  shown  in  DK83,  the  2DNH  is  initialized 

238  using  the  same  initial  conditions  and  mountain  profile  in  DK83  and  we  analyze  our  results 
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239 

240 

241 


242 

243 

244 

245 

246 

247 

248 

249 

250 

251 

252 

253 

254 

255 

256 

257 

258 


using  the  same  metrics  of  DK83. 

The  mountain  profile  is  given  by 

/?(x)  = - ^ - 2  (22) 

1+  - ^ 

V  y 

where  the  half-length  of  the  mountain  is  10  km,  the  height  is  1  m,  and  the 

prescribed  center  of  the  profile  is  0  km.  The  Initial  temperature  is  =  250  K  for  an 
isothermal  atmosphere  with  the  uniform  zonal  wind  u  =20  m/s.  In  the  isothermal  case,  the 

Brunt- Vaisala  frequency  =  g  ^  ~  A —  yields  the  potential  temperature  given  as 

dz  c 

P  0 


e 


(23) 


which  is  one  of  the  prognostic  variables  in  our  model.  The  domain  is  defined  as 
G  [-300,300]  X  [0,30]  km^.  The  bottom  boundary  uses  a  no-flux  boundary 

condition  while  the  lateral  and  top  boundaries  use  sponge  layers.  The  sponged  zone  is  10  km 
deep  from  the  top  and  50  km  wide  from  the  lateral  boundaries.  Over  the  sponge  layer  zone, 
the  prognostic  variables  are  relaxed  to  the  basic  initial  hydrostatic  state.  The  model  is 


integrated  for  a  nondimensional  time  of 


which  corresponds  to  8.33  hours  without 


diffusion  or  viscosity. 

Fig.  2  shows  the  numerical  and  analytic  solutions  at  steady  state  for  the  horizontal  and 
vertical  velocities,  which  agree  reasonably  well.  The  vertical  velocity  fields  match  very 
closely,  although  the  extrema  in  the  horizontal  velocity  field  are  underestimated  by  the 
numerical  model.  The  underestimated  extrema  in  the  horizontal  velocity  is  also  shown  in  both 
models  of  DK83  and  Giraldo  and  Restelli  (2008)  (hereafter,  GR08).  And  our  result  in  the 
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259  horizontal  velocity  is  in  good  agreement  with  DK83  and  GR08. 

260  Fig.  3  shows  the  normalized  momentum  flux  values  at  various  times  to  check  vertical 

261  transport  of  horizontal  momentum.  It  is  observed  that  the  flux  is  developing  well  and  the 

262  simulations  have  reached  steady-state  after  —  =  60  .  It  is  noted  that  the  mean  momentum 

a 

263  flux  at  that  time  is  97%  of  its  analytic  value.  It  agrees  well  with  DK83  as  well  as  GR08;  it  is 

264  important  to  point  out  that  the  Durran-Klemp  model  is  based  on  the  FD  method  in  both 

265  directions  while  the  Giraldo-Restelli  model  is  based  on  SEM  in  both  directions.  The 

266  mountain  test  shows  the  terrain-following  vertical  coordinate  is  well  suited  for  the 

267  combination  of  the  horizontal  SEM  and  vertical  EDM  for  spatial  discretization  even  though 

268  we  consider  a  small  mountain. 

269 

270  b.  2D  density  current  test 

271  In  order  to  verify  the  2DNH’s  feasibility  to  control  oscillations  with  numerical  viscosity 

272  and  evaluate  numerical  schemes  in  the  2DNH,  we  conduct  the  density  current  test  suggested 

273  by  Straka  et  al.  (1993).  The  density  current  test  is  initialized  using  a  cold  bubble  in  a  neutrally 

274  stratified  atmosphere.  When  the  bubble  touches  the  ground,  the  density  current  wave  starts  to 

275  spread  symmetrically  in  the  horizontal  direction  forming  Kelvin-Helmholtz  rotors.  Eollowing 

276  Straka  et  al.  (1993)  we  employ  a  dynamic  viscosity  of  v  =  75  m^s’'  to  obtain  converged 

277  numerical  solutions. 

278  Eor  an  initial  cold  bubble,  the  potential  temperature  perturbation  is  given  as 

Q 

279  6'  =  -^[_'\  +  cos{7rr)],  (24) 

with  the  center  of  the  bubble  at 


13 


281  (x^,zj  =  (0.3000)  m.  No-flux  boundary  conditions  are  used  for  all  boundaries.  The  model 

282  is  integrated  for  900  s  on  a  domain  -25600, 25600^  x  |^0, 6400^  m^. 

283  Fig.  4  shows  the  potential  temperature  perturbation  after  900s  for  400,  200,  100,  and 

284  50m  grid  spaeing  (Ax)  using  5th  order  basis  polynomials  per  element.  All  simulations  use 

285  Ai”  =  64  m  grid  spaeing  vertically.  As  expected,  the  higher  resolution  experiments  produee 

286  better  solutions  than  the  lower  resolution.  At  the  very  lowest  resolution  of  400  m,  only  two  of 

287  the  three  Kelvin-Helmholtz  rotors  are  generated  with  somewhat  eoarsened  frontal  surfaces.  In 

288  the  experiment  with  200  m  resolution,  the  three  rotors  appear  but  the  numerieal  solution  still 

289  suffers  from  coarsening  frontal  surfaees.  The  solutions  on  grids  finer  than  100  m  converge 

290  with  the  three  rotor  structures  adequately  simulated.  The  converged  solution  is  almost 

291  identical  to  other  published  solutions  (e.g.  Straka  et  al.  1993;  Skamaroek  and  Klemp  2008; 

292  GR08). 

293  In  order  to  show  the  effeet  of  higher  order  of  the  basis  polynomials,  we  show  the 

294  potential  temperature  perturbations  using  8th  order  basis  polynomials  per  element  with  the 

295  same  number  of  degrees  of  freedom  (DOF)  of  the  simulations  using  5th  order  basis 

296  polynomials  in  Fig.  5.  The  simulation  with  8th  order  basis  polynomials  on  the  very  lowest 

297  resolution  of  400  m  reprodueed  the  converged  solution  more  elosely  than  with  5th  order  basis 

298  polynomials.  Even  in  the  experiment  with  200  m  resolution,  the  coarsening  frontal  surfaees 

299  are  mitigated  and  the  solution  is  similar  with  the  eonverged  solution  with  three  rotors. 

300  Fig.  6  shows  the  profiles  of  the  potential  temperature  perturbation  at  the  height  of  1200 

301  m.  The  results  from  the  highest  grid  resolution  of  the  simulations  using  5th  and  8th  order 

302  basis  polynomials  are  indistinguishable  and  well  converged  (Fig.  6a).  Three  minima 

303  eorresponding  to  the  three  rotors  agree  well  with  other  published  solutions.  In  addition  to  the 

304  profiles,  the  front  loeation  (-1K  of  potential  temperature  perturbation  at  the  surface),  and  the 
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extrema  of  the  pressure  perturbation  and  potential  temperature  perturbation  agree  well  with 
eaeh  other  (Table  1),  of  whieh  the  numbers  are  comparable  to  those  of  GR08.  In  the 
numerical  results  from  the  different  grid  resolutions  simulated  by  using  5th  order  basis 
polynomials,  the  potential  temperature  profile  at  the  coarsest  resolution  of  400  m  grid  shows 
significant  fluctuations  (Fig.  6b).  That  of  8th  order  polynomials,  however,  tends  to  be 
relieved  from  the  deviation  from  the  converged  solution  (Fig.  6c).  The  above  results  suggest 
that  the  numerical  solution  can  be  converged  more  rapidly  by  using  higher  order  of  basis 
polynomial.  Furthermore,  the  results  in  this  paper  show  that  an  adequate  convergence  can  be 
reached  at  grid  resolutions  finer  than  200  m. 


c.  Inertia-gravity  wave  test 

This  test  examines  the  evolution  of  a  potential  temperature  perturbation  6'  in  a 
constant  mean  flow  with  a  stratified  atmosphere.  This  initial  perturbation  diverges  to  the  left 
and  right  symmetrically  in  a  channel  with  periodic  lateral  boundary  conditions.  The  inertia- 
gravity  wave  test  introduced  by  Skamarock  and  Klemp  (1994)  (hereafter,  SK94)  serves  as  a 
tool  to  investigate  the  accuracy  for  NH  dynamics.  Also  we  use  this  experiment  to  check  the 
consistency  of  the  results  with  various  resolutions.  The  parameters  for  the  test  are  the  same  as 
those  of  SK94.  The  initial  state  is  a  constant  Brunt- Vaisala  frequency  of  N  =  0.01/s  with 
surface  potential  temperature  of  =  300  K  and  a  uniform  zonal  wind  U  -  20  m/s.  In  order 
to  trigger  the  wave,  the  initial  potential  temperature  perturbation  9'  is  overlaid  to  the  above 
initial  state  and  is  given  as 


Sin 


e'(x.z)  =  e^ 


1  + 


V  y 


(25) 
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where  ^=0.01K,  z  =10  km,  x  =100  km.  The  domain  is  defined  as 

328  ^  [0,300]  X  [0,1 0]  km^.  We  use  periodic  lateral  boundary  conditions  and  a  no-flux 

329  boundary  conditions  for  both  the  bottom  and  top  boundaries.  The  simulation  is  performed  for 

330  3000s  with  no  viscosity. 

331  Fig.  7  shows  the  solution  6'  at  the  initial  time  and  time  3000  s  with  a  horizontal 

332  resolution  Ax  =  250  m  and  a  vertical  resolution  Az  =  250  m.  The  figure  uses  the  same 

333  contouring  interval  as  in  SK94  and  Giraldo  and  Restelli  (2008)  for  comparison.  The  results 

334  are  produced  with  8th  order  polynomials  per  element.  We  have  conducted  the  2DNH  model 

335  with  various  basis  polynomial  orders  at  the  same  resolution,  where  the  simulated  results  are 

336  found  to  be  very  comparable.  SK94  give  an  analytic  solution  for  the  case  of  the  Boussinesq 

337  equations,  but  it  is  only  valid  for  the  Boussinesq  equations  while  we  use  the  fully 

338  compressible  equations  in  our  model.  Using  the  analytic  solution  only  for  qualitative 

339  comparisons,  we  find  that  the  extrema  of  our  results  are  comparable  to  the  analytic  values.  In 

340  comparison  with  the  results  of  Giraldo  and  Restelli  (2008)  in  which  the  fully  compressible 

341  equations  are  also  used,  our  results  look  very  similar.  Fig.  8  shows  the  profiles  along  5000  m 

342  for  various  horizontal  resolutions.  All  models  show  consistently  identical  solutions  with  the 

343  symmetric  distribution  about  the  midpoint  ( x  =  1  60  km)  which  is  the  location  to  which  the 

344  initial  perturbation  moved  by  the  horizontal  flow  of  20  m/s  after  3000  s.  Even  at  coarser 

345  resolution  experiments,  it  does  not  exhibit  phase  errors  although  the  maxima  and  minima  near 

346  the  midpoint  (x  =  160km)  are  slightly  damped.  Table  2  shows  the  extrema  of  vertical 

347  velocities  and  potential  temperature  perturbations  for  the  results  of  various  horizontal 

348  resolutions  after  3000  s.  It  is  noted  that  all  experiments  give  almost  the  same  values  for 

349  potential  temperature  perturbation  where  these  values  in  the  range 

350  0'  e  1^-1 .52  X 1 0“^,2.83  x  1 0“^  j  are  comparable  to  other  studies  (e.g.,  GR08  and  Li  et  al. 
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351  2013).  GR08  give  the  ranges  of  6'  g  -1 .51  x  1  0  ^,2.78  x  1 0  ^  from  the  model  based  on 


352  the  spectral  element  and  discontinuous  Galerkin  method.  Also  Li  et  al.  (2013)  show 

353  G  1^-1. 53  xlO“^,  2.80  xlO”^J  using  the  high-order  conservative  finite  volume  model 

354  which  are  similar  to  our  results. 

355 

356  d.  Rising  thermal  bubble  test 

357  We  also  conduct  the  rising  thermal  bubble  test  to  verify  the  consistency  of  the  scheme  in 

358  the  model  to  simulate  thermodynamic  motion  (Wicker  and  Skamarock  1998).  This  test 

359  considers  the  time  evolution  of  warm  air  in  a  constant  potential  temperature  environment  for 

360  an  atmosphere  at  rest  atmosphere.  The  air  that  is  warmer  than  the  ambient  air  rises  due  to 

361  buoyant  forcing  which  then  deforms  due  to  the  shearing  motion  caused  by  gradients  of  the 

362  velocity  field  and  eventually  shapes  the  thermal  bubble  into  a  mushroom  cloud.  Because  the 

363  test  case  has  no  analytic  solution,  the  simulation  results  are  evaluated  qualitatively. 

364  The  initial  conditions  we  use  follow  those  of  GR08  in  which  the  domain  for  the  case  is 

365  defined  as  (x,z)g[0,1]  km^.  We  consider  no-fiux  boundary  conditions  for  all  four 


366  boundaries.  The  domain  is  initialized  for  a  neutral  atmosphere  at  rest  with  6^  =  300  K  in 

367  hydrostatic  balance.  A  potential  temperature  perturbation  to  drive  the  motion  is  given  as 


368 


for  r  >  r 

c 

for  r  <r 


(26) 


369  where  0^  =0.5K,  r  =  J{x  with  ^  =  (SOO,  350^m.  The  model 


370  was  run  for  a  time  of  700  seconds.  It  should  be  noted  that  an  explicit  second-order  diffusion 

371  on  coordinate  surfaces  is  used  with  a  viscosity  coefficient  of  v  =  1  m^s"'  for  all  simulations 
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372  of  this  test.  The  numerical  diffusion  is  applied  for  momentum  and  potential  temperature  along 

373  the  horizontal  and  vertical  directions  so  that  it  eliminates  the  erroneous  oscillations  at  the 

374  small  scale  -  while  this  amount  of  diffusion  might  seem  excessive,  it  has  been  chosen 

375  because  it  allows  the  model  to  remain  stable  even  after  the  bubble  hits  the  top  boundary. 

376  Fig.  9  shows  the  potential  temperature  perturbation,  horizontal  wind,  and  vertical  wind 

377  fields  for  the  simulations  of  two  resolutions  of  20  m  and  5  m  horizontal  and  vertical  grid 

378  spacings  ( Ax  and  Ai~ ),  respectively,  employing  5th  order  basis  polynomials.  In  both 

379  simulations,  the  line  structures  in  the  numerical  solutions  are  well  depicted  with  a  perfectly 

380  symmetric  distribution  at  the  midpoint  and  sharp  discontinuities  of  the  fields  along  boundary 

381  lines  of  the  bubble.  At  lower  resolution,  however,  degradations  in  the  solution  are  visible  in 

382  the  potential  temperature  perturbation  and  vertical  wind  which  are  illustrated  by  fluctuations 

383  in  the  values  as  well  as  the  concaving  contours  at  the  top  of  the  bubble.  It  is  noted  that  while 

384  the  numerical  solution  of  the  model  using  the  spatially  centered  FDM  of  Wicker  and 

385  Skamarock  (1998)  shows  spurious  oscillations  in  the  potential  temperature  field,  the 

386  simulations  here  of  2DNH  using  the  horizontally  SEM  and  vertically  FDM  is  devoid  of  these 

387  oscillations. 

388  We  also  show  the  vertical  profiles  of  potential  perturbation  at  x  =  500  m  after  700  s  for 

389  various  resolutions  in  Fig.  10.  Simulations  were  run  with  various  resolutions  of  5,  10,  and  20 

390  m,  where  the  resolutions  given  are  defined  for  both  the  horizontal  and  vertical  directions.  The 

391  results  of  10  m  and  5  m  resolutions  are  almost  identical  to  each  other.  The  result  of  the  lowest 

392  resolution  of  20  m,  however,  shows  a  somewhat  unresolved  solution,  in  which  the  maximum 

393  value  is  underestimated  and  the  phase  shift  is  depicted.  The  time  series  for  maximum 

394  potential  temperature  perturbation  and  maximum  vertical  velocity  are  shown  in  Fig.  1 1.  In  all 

395  simulations,  the  maximum  vertical  velocity  increases  as  the  maximum  theta  perturbation 

396  decreases.  This  shows  that  the  thermal  energy  of  the  theta  perturbation  leads  to  the 
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397  acceleration  of  the  vertieal  veloeity.  This  result  agrees  well  with  the  study  of  Ahmad  and 

398  Lindeman  (2007). 

399 

400  6.  Summary  and  Conclusions 

401  The  non-hydrostatie  eompressible  Euler  equations  for  a  dry  atmosphere  are  solved  in  a 

402  simplified  2D  sliee  (X-Z)  framework  by  using  the  speetral  element  discretization  (SEM)  in 

403  the  horizontal  and  the  third-order  finite  differenee  seheme  for  the  vertieal  diseretization.  The 

404  form  of  the  Euler  equations  used  here  are  the  same  as  those  used  in  the  Weather  Researeh  and 

405  Eoreeasting  (WRE)  model.  We  employ  a  hybrid  sigma-pressure  vertieal  eoordinate  whieh  ean 

406  be  modified  exactly  to  the  sigma-pressure  eoordinate  at  the  level  of  the  actual  coding 

407  implementation. 

408  Eor  the  spatial  discretization,  the  spatial  operators  are  separated  into  their  horizontal  and 

409  vertieal  eomponents.  In  the  horizontal  eomponents,  the  operators  are  discretized  using  the 

410  SEM  in  whieh  high-order  representations  are  construeted  through  the  GEE  grid  points  by 

411  Eagrange  interpolations  in  the  elements.  Eising  GEE  points  for  both  interpolation  and 

412  integration  results  in  a  diagonal  mass  matrix,  whieh  means  that  the  inversion  of  the  mass 

413  matrix  is  trivial.  In  the  vertieal  eomponents,  the  operators  are  diseretized  using  the  third-order 

414  upwind  biased  finite  differenee  scheme  for  the  vertieal  fluxes  and  eentered  differenees  for  the 

415  vertieal  derivatives.  The  time  diseretization  relies  on  the  time-split  third-order  Runge-Kutta 

416  teehnique. 

417  We  have  presented  idealized  standard  benehmark  tests  for  large-seale  flows  (e.g.,  linear 

418  hydrostatie  mountain  wave)  and  for  nonhydrostatie-seale  flows  (e.g.,  inertia-gravity  wave, 

419  rising  thermal  bubble,  and  density  eurrents).  The  numerieal  results  show  that  the  present 

420  dynamieal  eore  is  able  to  produee  solutions  of  good  quality  eomparable  to  other  published 
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421  solutions.  These  tests  effeetively  reveal  that  the  eombined  spatial  diseretization  method  of  the 

422  spectral  element  and  finite  difference  method  in  the  horizontal  and  vertical  directions, 

423  respectively,  offers  a  viable  method  for  the  development  of  a  NH  dynamical  core.  Further 

424  research  will  be  continued  to  couple  the  present  core  with  the  existing  physics  packages 

425  together  and  extend  the  2D  slice  framework  to  develop  a  3D  dynamical  core  for  the  global 

426  atmosphere  where  the  cubed-sphere  grid  will  be  used  for  the  spherical  geometry. 

427 
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533  Table  Captions 

534  Table  1.  Comparison  between  the  5th  and  8th  order  polynomials  per  elements  for  the 

535  density  eurrent.  The  simulations  is  eondueted  with  Ax  =50  m  and  Az~  =  50  m  resolution. 

536 

537  Table  2.  Comparison  of  the  numerieal  results  for  various  horizontal  resolutions.  All 

538  simulations  use  the  8th  order  polynomials  per  elements  and  vertieal  resolution  of 

539  AT  =  250  m. 

540 
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541  Table  1.  Comparison  between  the  5th  and  8th  order  polynomials  per  elements  for  the 

542  density  current.  The  simulations  is  condueted  with  Ax  =50  m  and  Az~  =  50  m  resolution. 


543 

544 

545 

546 

547 


548 


Order  of 
polynomials 

Front 

location(km) 

P '  (Pa) 

P '  (Pa) 

~  min^  ' 

e'  (K) 

max  ^  ' 

e' .  (K) 

mm  ^  ' 

5  th 

14.77 

630.62 

-452.79 

0.08 

-8.87 

8th 

14.74 

626.91 

-456.84 

0.08 

-8.94 

Table  2.  Comparison  of  the  numerical  results  for  various  horizontal  resolutions  for 
inertia-gravity  wave.  All  simulations  use  the  8th  order  polynomials  per  elements  and  vertical 


resolution  of  Az  =  250m. 


Resolution(m) 

W  (m/s) 

1/1/  .  (m/s) 

mm  ^  ' 

0'  (K) 

max  ^  ' 

0'  (K) 

mm  ^  ' 

Ax  =125 

2.85X10'^ 

-2.89X10'^ 

2.83X10'^ 

-1.52X10'^ 

Ax  =250 

2.80X10'^ 

-2.82X10'^ 

2.83X10'^ 

-1.52X10'^ 

Ax  =500 

2.73X10'^ 

-2.73X10'^ 

2.83X10'^ 

-1.52X10'^ 

Ax  =750 

2.72X10'^ 

-2.70X10'^ 

2.83X10'^ 

-1.52X10'^ 

Ax  =1250 

2.68X10'^ 

-2.62X10'^ 

2.82X10'^ 

-1.52X10'^ 

549 


27 


550  Figure  Captions 


551  FIG.  1.  The  grid  points  of  columns  within  an  element  having  four  GLL  points.  The 

552  hybrid  sigma  coordinate  are  illustrated  and  the  close  (open)  circles  on  the  solid  (dashed)  line 

553  indicate  the  location  of  the  variables  at  layer  mid-points  (interfaces). 

554 

555  FIG.  2.  Steady-state  flow  of  (left)  horizontal  velocity  (m/s)  and  (right)  vertical  velocity 

556  (m/s)  over  1  m  high  mountain  at  nondimensional  time  —  =  60  with  a  grid  resolution  of 

a 

557  Ax  =  2  km  using  5th  order  basis  polynomials  per  element  and  Az~  =  375  m.  The 

558  numerical  solution  is  represented  by  solid  lines  and  the  analytic  solution  by  dashed  lines. 

559 

560  FIG.  3.  Vertical  flux  of  horizontal  momentum,  normalized  by  its  analytic  value  at 

561  several  non-dimensional  times  — .  Here  M  and  Mh  are  the  momentum  flux  of  the  numerical 

a 

562  and  analytic  solution. 

563 

564  FIG.  4.  Potential  temperature  perturbation  after  900  s  using  (a)  Ax  =  400  m,  (b) 

565  Ax  =  200  m,  (c)  Ax~  =100m,  and  (d)  Ax  =50m  grid  spacing  with  5th  order  basis 

566  polynomials  per  element.  All  simulations  use  Az~  =  64  m  grid  spacing. 

567 

568  FIG.  5.  As  in  Fig.  4,  but  with  8th  order  basis  polynomials  per  element. 

569 

570  FIG.  6.  Profiles  of  potential  temperature  perturbation  after  900  s  along  1200  m  height: 

571  (a)  high-resolution  simulations  using  5th  (thin  solid  line)  and  8th  (thick  solid  line)  order  basis 

572  function,  (b)  simulations  using  5th  order  basis  polynomials,  and  (c)  simulations  using  8th 
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573  order  basis  polynomials  per  element.  Total  number  of  the  degrees  of  freedom  is  the  same  in 

574  both  5th  and  8th  order  experiments.  All  simulations  use  Ai”  =  64  m  grid  spaeing. 

575 

576  FIG.  7.  Potential  temperature  perturbation  at  the  initial  time  (left)  and  time  3000s  (right) 

577  for  Ax  =  250  m  using  8*  order  basis  polynomials  per  element  and  Az~  =  250  m. 

578 

579  FIG.  8.  The  profiles  of  potential  temperature  perturbation  along  5000  m  height  for 

580  Ax  =125m  (thick  solid  line),  Ax  =500m  (thin  dashed  line)  and  Ax"  =  1250  m  (thin 

581  solid  line)  using  8th  order  basis  polynomials  per  element.  All  models  use  Ai”  =  250  m. 

582 

583  FIG.  9.  Plots  of  (a,b)  potential  temperature  perturbation  (K),  (c,d)  horizontal  wind  (m/s), 

584  and  (e,f)  vertical  wind  (m/s)  for  the  rising  thermal  bubble  test  after  700s  with  (left) 

585  Ax”  =  20  m  and  (right)  Ax”  =  5  m  resolution.  All  simulations  use  5 order  basis  polynomials 

586  per  element  and  Ai”  =  10  m  grid  spacing.  All  negative  values  are  denoted  by  dashed  lines 

587  and  positive  values  by  solid  lines. 

588 

589  FIG.  10.  Vertical  profiles  of  the  potential  temperature  perturbation  at  x  =  500  m  after 

590  700  s  for  various  resolutions:  Ax^,  Ai”  =  20  m  (thin  solid  line).  Ax”,  Ai”  =  10m  (thin  dashed 

591  line),  and  Ax^,  Ai”  =  5  m  (thick  solid  line). 

592 

593  FIG.  II.  (top)  Domain  maximum  potential  temperature  perturbation  and  (bottom) 

594  vertical  wind  for  the  rising  thermal  bubble  test.  All  simulations  use  the  5th  order  basis 

595  polynomials  per  element,  and  the  vertical  resolutions  are  the  same  as  the  horizontal 

596  resolutions. 

597 
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FIG.  1.  The  grid  points  of  columns  within  an  element  having  four  GLL  points.  The 
hybrid  sigma  coordinate  are  illustrated  and  the  close  (open)  circles  on  the  solid  (dashed)  line 
indicate  the  location  of  the  variables  at  layer  mid-points  (interfaces). 
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FIG.  2.  Steady-state  flow  of  (left)  horizontal  velocity  (m/s)  and  (right)  vertical  velocity 

(m/s)  over  1  m  high  mountain  at  nondimensional  time  —  =  60  with  a  grid  resolution  of 

a 

Ax  =  2  km  using  5th  order  basis  polynomials  per  element  and  AT  =  375  m.  The 
numerical  solution  is  represented  by  solid  lines  and  the  analytic  solution  by  dashed  lines. 
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608 

609 

610 

611 

612 


Vertical  flux  of  horizontal  momentum 


FIG.  3.  Vertical  flux  of  horizontal  momentum,  normalized  by  its  analytic  value  at 

several  non-dimensional  times  — .  Here  M  and  Mh  are  the  momentum  flux  of  the  numerical 

a 

and  analytic  solutions. 
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613  X(km) 

614  FIG.  4.  Potential  temperature  perturbation  after  900  s  using  (a)  Ax  =  400  m,  (b) 

615  Ax  =  200  m,  (e)  Ax  =100m,  and  (d)  Ax  =  50  m  grid  spacing  with  5th  order  basis 

616  polynomials  per  element  for  the  density  current.  All  simulations  use  Az~  =  64m  grid 

617  spacing. 
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620  FIG.  5.  As  in  Fig.  4,  but  with  8th  order  basis  polynomials  per  element. 
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FIG.  6.  Profiles  of  potential  temperature  perturbation  after  900  s  along  1200  m  height: 

(a)  high-resolution  simulations  using  5th  (thin  solid  line)  and  8th  (thick  solid  line)  order  basis 

function,  (b)  simulations  using  5th  order  basis  polynomials,  and  (c)  simulations  using  8th 

order  basis  polynomials  per  element.  The  total  number  of  the  degrees  of  freedom  is  the  same 

in  both  5th  and  8th  order  experiments.  All  simulations  use  Ai”  =  64  m  grid  spacing. 
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629  X  (km)  X  (km) 

630  FIG.  7.  Potential  temperature  perturbation  at  the  initial  time  (left)  and  time  3000s  (right) 

631  for  Ax  =  250  m  using  8*’’  order  basis  polynomials  per  element  and  AT  =  250  m  for  the 

632  inertia-gravity  wave. 
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FIG.  8.  The  profiles  of  potential  temperature  perturbation  along  5000  m  height  for 
Ax  =125m  (thick  solid  line),  Ax  =  500m  (thin  dashed  line)  and  Ax  =  1250m  (thin 
solid  line)  using  8th  order  basis  polynomials  per  element  for  the  inertia-gravity  wave.  All 
models  use  Az~  =  250m. 
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641  FIG.  9.  Plots  of  (a,b)  potential  temperature  perturbation  (K),  (o,d)  horizontal  wind  (m/s), 

642  and  (e,f)  vertieal  wind  (m/s)  for  the  rising  thermal  bubble  test  after  700s  with  (left) 

643  Ax,  AT  =  20 m  and  (right)  Ax,  Az  =  5  m  resolution  for  the  rising  thermal  bubble  test.  All 

644  simulations  use  5'*'  order  basis  polynomials  per  element.  All  negative  values  are  denoted  by 


645  dashed  lines  and  positive  values  by  solid  lines. 
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FIG.  10.  Vertical  profiles  of  the  potential  temperature  perturbation  for  the  rising  thermal 
bubble  test  at  x  =  500  m  after  700  s  for  various  resolutions;  Ax,  AT  =  20 m  (thin  solid  line), 
Ax,  Ai”  =  10m  (thin  dashed  line),  and  Ax,  Ai~  =  5  m  (thick  solid  line). 
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654  FIG.  11.  (top)  Domain  maximum  potential  temperature  perturbation  and  (bottom) 

655  vertieal  wind  for  the  rising  thermal  bubble  test.  All  simulations  use  the  5th  order  basis 

656  polynomials  per  element,  and  the  vertieal  resolutions  are  the  same  as  the  horizontal 

657  resolutions. 
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